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We show that with a small modification, the formulation of the Einstein equations of Uggla et al, 
which uses tetrad variables normalised by the expansion, is a mixed symmetric hyperbolic/parabolic 
system. Well-posedness of the Cauchy problem follows from a standard theorem. 
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I. INTRODUCTION 

There has been much interest recently in writing the 
Einstein field equations in different ways and analyzing 
the properties of the resulting equations. Much of the 
motivation for this comes from the search for a form of 
the equations suitable for numerical simulations. In all 
these different forms, the Einstein field equations form a 
constrained system. Since numerical evolutions can only 
approximate the constraints, an analysis of the equations 
must look at the larger solution space where the con- 
straints arc not satisfied. In the larger solution space, 
different formulations of the field equations can have very 
different properties, in terms both of well-posedness and 
of stability. 

One unusual form of the field equations is due to Ug- 
gla et al^ who use a set of scale invariant variables 
based on the tetrad formalism. The motivation of 
was not to find a system suitable for numerical simula- 
tions; but rather to find a system suitable for the anal- 
ysis of the properties of singularities, and in particular 
of the asymptotic properties of the metric as the singu- 
larity is approached. For this purpose, well-posedness is 
not needed: it is well known that there are forms of the 
Einstein field equations for which the initial value formu- 
lation is well-posed. One can then take solutions of the 
field equations (in one of these forms) and ask what the 
behavior of the variables of 1\ is in these solutions as the 
singularity is approached. 

Nonetheless, one can examine the equations of lll| to 
see whether they have the properties that one would need 
{e.g. well-posedness) for the system to be suitable for 
numerical simulations. Here, the situation does not look 
hopeful. A cursory examination of the principal part of 
the equations of 0] shows that it is not of any standard 
type (hyperbolic or parabolic) that would lead one to 
expect well-posedness. The equations of were refor- 
mulated in '2?| and then used for a numerical simulation 
of the approach to the singularity. This reformulation es- 
sentially involved using an evolution equation for a quan- 
tity that had previously been given by a constraint. In 
the form used in 2\ the principal part of the equations 
has the form of a combination of advection and diffusion 
equations. Therefore, one might hope that the system 



was well-posed. However, no analysis of well-posedness 
was done in 0. 

Here, we perform such an analysis. We show that the 
system of |3] is not well-posed as written, but that it 
can be made so by a simple modification. The modifica- 
tion consists of adding a multiple of one of the constraint 
equations to one of the evolution eqautions. 

In section II we introduce the systems of and bi- 
section III contains the analysis of the well-posedness of 
the main system, and Section IV analyses the implied 
constraint evolution system. Conclusions are given in 
section V. 



II. EQUATIONS 

In P] spacetime is described in terms of a set of coordi- 
nates (t, a;') and a set of orthonormal vectors (the tetrad) 
(gq, Ba) where both the spatial coordinate index i and the 
spatial triad index a go from 1 to 3. The timelike vector 
Go of the tetrad is chosen to be hypersurface orthogonal 
with the relation between tetrad and coordinates chosen 
to be of the form Gq — N~^dt and g„ — Cadi. The 
commutators of the tetrad components are decomposed 
as follows: 
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(1) 
(2) 



where n"-^ is symmetric, cr"'^ is symmetric and trace free, 
and Cap-y is totally antisymmetric with £123 = 1. 

Scale invariant variables are defined by dividing ap- 
propriate quantities by the expansion H . In particular, 
scale invariant frame derivatives are given by {da, da} = 
{gq, Ba} / H while scale invariant derivatives of the expan- 
sion are given by g -f 1 = — 9o Ini? and = —da IniJ. 
The lapse N is chosen to be equal to . This is equiv- 
alent to the statement that the surfaces of constant time 
are chosen to evolve according to inverse mean curva- 
ture flow. This condition yields do — dt and iia — Hra- 
Other scale invariant variables are given by 

{Ea\^ap,A'',Na(i} = {e„\ (7,^, fl", n„^}/i/ (3) 

In Q the dynamical system is given by the scale invari- 
ant variables of equation along with . The quantity 
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is not a dynamical variable as it can be chosen arbi- 
trarily by choosing a rule for the propagation of the spa- 
tial triad . In particular, if the spatial triad is chosen 
to be Fermi propagated, then fl" = 0. 

The quantity q is also not a dynamical variable in (ij 
as it is given by 

q - iE^/'E,/, - ia„r" -t- f ^„r" (4) 

where equation Q| follows from the vacuum Einstein 
equation. Einstein's equation then becomes a set of evo- 
lution and constraint equations for the dynamical vari- 
ables. The evolution equations are all first order in time, 
and with the exception of the evolution equation for 
they are all first order in space. The evolution equation 
for Ta can be written as 

dtro. = a„g + (qSj - + e^^^K^^p (5) 

where i?" = il'^/H. However, since q is not one of the 
dynamical variables, the quantities q and daq must be 
expressed by solving equation Q for q. The result is 

+ (iE^^E^^ - + lApr^) r„ (6) 

The first term on the right hand side of equation © is 
the one with the highest number of derivatives. It is 
second order, but not positive definite, so the equation 
is not parabolic. It does not seem to be of a type for 
which results about well-posedness are known and does 
not seem to be the sort of equation for which a numerical 
evolution is likely to be stable. 

To address these issues, a slight variation of the system 
of was used in Here, the quantity q was added to 
the list of dynamical variables. The evolution equation 
for Ta is then simply equation ((SJ and is first order in 
space. Equation |@J then becomes a constraint equation. 
Since q is now a dynamical variable, it has an evolution 
equation which turns out to be second order in space. 
The full set of evolution equations for the system of ,2j is 

dtE^' = FjEp' (7) 
dtvc. = Fjrp+d^q (8) 
a^A" = F'^pA^ + \dpT.°'^ (9) 

+ a<"r'3> - a<"A''> + 2r<"A^> (10) 
+ e'^<"{d^-2A^)N^>s (11) 
dtN""^ = qN''^ + 2Y.''"sNP^^ ~ e^^^^'d^Y.^^ (12) 
dtq = [2{q - 2) + i (2A" - r") d„ - \d"d^] q 

- |a„r" -f |A"r„ + |r/3d„E"'3 

- 2E"'5VK„^ (13) 

Here angle brackets denote the symmetric trace- free part, 
and Fap and Waf} are given by 

Fap = qSaP - Ya0 (14) 



- ldo,rp-\e'^^^{d-,-2A^)Nps (15) 

The gauge choice fi" = is made, which means that the 
spatial triad is Fermi propagated. Note that the highest 
order term in equation (|13|l is —{l/'i)8Pdaq which resem- 
bles the backward diffusion equation. Thus, we might 
expect this equation to be well behaved for evolution in 
the negative time direction, which is the direction appro- 
priate for the approach to the singularity. 

In addition to the evolution equations, the variables 
satisfy constraint equations as follows: 

= {Ccorafap = 2(5[q - r[„ - A[a)Ei3] 

- e^psN^^E^' (16) 
- Cg = 1 + i(2a„ - 2r, - 3A,)A" - ^N^pN'^'^ 

+ UN-cf ~ l^cpY'^' (17) 
= {CcT =dpT.°'^ + 2r°' -Y."prP -iApY.'^'^ 

- e^^-'Npsi:^' (18) 
= Cq = q- iE^'^E^/s + ia„r" - |A„r" (19) 
= (Cj)" EE {dp - rp^N'^" + e"^^A^) 

- 2ApN'''^ (20) 
= {CwT = [e^f'^dp - Ap) - 7V"^]r^ (21) 

III. ANALYSIS OF WELL-POSEDNESS 

We now analyze the evolution equations of Q (equa- 
tions (|7I13|) to see whether they are well-posed. We shall 
show that with a small modification, the system can be 
brought into a form that is mixed parabolic/hyperbolic. 
We then use a textbook theorem to show that the result- 
ing initial value problem (without boundaries) is well- 
posed. 

The relevant theorems are theorems 4.6.2 and 4.9.1 of 
. Here the first theorem is for linear systems, while the 
second theorem applies the results of the first theorem to 
nonlinear systems by considering them as a sequence of 
linear systems. For simplicity, we will present the anal- 
ysis of a single linear system in the sequence. That is 
we assume a guess for the dynamical variables, use that 
guess in the coefficients of the evolution equations H7I13() 
and consider the resulting linear problem. 

Theorem 4.6.2 of asserts the following: Assume we 
have a vector of variables u and another vector of vari- 
ables V, which obey a linear system of evolution equations 
of the form 

dtu = Diiu + Di2V, (22) 
dtv = D21U + D22V. (23) 

Here the D are linear spatial derivative operators whose 
coefficients can depend on t and x^. Du is a first-order 
derivative operator such that dtu = Duu is symmet- 
ric hyperbolic. D22 is a second-order derivative opera- 
tor such that dtV = D22V is parabolic. D12 and D21 
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are arbitrary first-order derivative operators. Then the 
coupled problem is said to be mixed symmetric hyper- 
bolic/parabolic, and is well-posed. The theorem assumes 
periodic boundary conditions, which is what we are in- 
terested in in applications to cosmology. 

Garfinkle's system is first order in time and first order 
in space, except for the second derivative in dtq = . . . — 
1/3 dad^q. We therefore try to apply the above theorem 
to 

v = q, u=(i?„',r",s",S"^iV"''), (24) 

Here to simplify the expressions for the differential op- 
erators we have introduced the notation s" = A" — r". 
Well-posedness depends only on the principal part (here, 
terms containing first derivatives) and in examining it we 
can assume that the coefficients are constant. In partic- 
ular, we can treat the frame derivatives 

a„^i?.'A (25) 

as if they were partial derivatives, because their commu- 
tator is a lower-order term. In the following, we use the 
symbol da to stress this. The principal part of D22 is 

dtq^-^d'^daq. (26) 

Here ~ stands for "equal up to lower order terms" . It is 
clear that dtv — D22V is parabolic (going backwards in t). 
Furthermore Du, D12 and D21 are first order. We only 
need to check that dtu — Duu is symmetric hyperbolic. 

As in any evolution system subject to constraints, we 
can change the level of hyperbolicity of the evolution 
equations for u by adding constraints to the right-hand 
side of the evolution equations. This does not affect so- 
lutions of the evolution equations which also obey the 
constraints, but changes the evolution equations off the 
constraint surface. The evolution equations are first or- 
der and so are the constraints. Therefore we can add any 
undifferentiated constraint to the right-hand side that 
has the correct number of indices and symmetries. In 
the current system, we can add arbitrary multiples of 
(Cc)", {Cj)" and (Cw)" to the evolution equations for 
r" and s" . We can also add 5"^ times arbitrary multiples 
of Cg and Cq to the evolution equation for N"^ . There 
are no other possibilities, and so addition of constraints 
is determined by eight free parameters. 

The complete calculation is straightforward but te- 
dious. For our purposes it will be sufficient to add d(Cc)" 
to dt-s", where d is a constant parameter, d = corre- 
sponds to the system of equations of . The principal 
part of Dii, making the symmetrizations explicit, is 

dtEa' ~ 0, (27) 
dtr" ~ 0, (28) 

dts" ^ (d+^\dpJ:"f', (29) 



+ ^ie^^"d^NPs + e^^^d^N^s), (30) 
dtN^'f^ ~ -^{e'-'^'^d^Y.'^s + e''^^d^T,°'s)- (31) 

In order to show symmetric hyperbolicity, we have to 
do two separate things. First we have to find a set of 
characteristic variables {U} that is complete in the sense 
of spanning {u}. Let n" be a unit vector with respect 
to the metric dap- Then a characteristic variable U with 
speed —A in the direction n" is a linear combination of 
u with the property that 

dtU = Xn°'daU + derivatives normal to n" 

+lower order terms. (32) 

If all speeds A are real, and the U are complete, the 
system is said to be strongly hyperbolic. We introduce a 
notation where an index n denotes contraction with n" or 
Ua, and where uppercase Latin indices denote projection 
with the projector 

9a/3 = 6ap - riaUp. (33) 

In this notation, the principal part of the evolution equa- 
tions is 

dtu = P'^daU = (P„9„ + P^dA)u. (34) 

The system is strongly hyperbolic if and only if P„ has 
only real eigenvalues and is diagonalizable, and if the 
eigenvalues and the diagonalizing matrix depend contin- 
uously on the direction n" . 

Secondly, we have to find an energy density e that is 
quadratic in {u}, conserved in the sense that there exists 
a flux so that 

dt€ = daF"' + lower order terms (35) 

and that is positive definite in the sense that it is non- 
negative, and zero if and only if m — 0. The system is said 
to be symmetric hyperbolic if and only if it is strongly 
hyperbolic and if it admits a positive definite conserved 
energy. 

In order to show strong hyperbolicity, we bring the ma- 
trix Pn into block-diagonal form by considering suitable 
linear combinations of the u with respect to n" . Any vari- 
able whose time derivative is zero and whose derivative in 
the n" direction does not appear in the evolution equa- 
tions is a characteristic variable with zero speed (called 
a zero mode). For the purposes of this analysis, each 
zero mode can be considered a separate 1 by 1 block of 
Pn that automatically satisfies the conditions needed for 
strong hyperbolicity. From equations H27I31|I it immedi- 
ately follows that Ea^ and r" are zero modes. 

We define s„ as the contraction of Sa with n", sa as 
its projection with qap, Nqq as N^^qap, and similarly for 
other variables. We similarly project the evolution equa- 
tions. Using these projections we find that 7V„„ and Nqq 
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are zero modes and that for the remaining variables P„ is 
block-diagonal in three blocks, corresponding to variables 
which are, respectively, scalars, vectors, and symmetric 
tracefree tensors in the space normal to n". 
The scalar block is given by 



5ti 



d 

-2/3 



-1/2 




(36) 



Note that because Sq^ is tracefree (with respect to Sap), 
we have E„„ -I- Sg, — 0. Therefore we do not need to 
consider Eg, as a separate variable. P.^ has eigenvalues 
±a/— (1 -I- 2c?)/3, and it is diagonalizable for d ^ —1/2. 
The vector block is given by 



(37) 













J ^ n 1 










Here we have defined the shorthand 

^An = enBAN^"^. 



(38) 



P,^ has eigenvalues (0, ±-y/— (i/2), and is diagonalizable 
for d ^ 0. 

The symmetric tracefree tensor block is given by 



(AB) 
^AB 



P' = 



1 

1 



We have defined the shorthands 



Eab = '^AB — -^QAB^qq, 
Nab = enCANe^. 



(39) 

(40) 
(41) 



Note that it is consistent to remove the g-trace from E as 
Yiqq — — E„„ appears in the scalar block. Note also that 
although the 5-trace of E^^ vanishes, this does not mean 
that the q-trace of the projection Y,ab vanishes. Finally, 
note that Nab has vanishing q-trace by definition, but 
is not symmetric. What appears in the tensor block is 
its symmetrisation N(^ab) ■ This is consistent because the 
antisymmetric part is N^^ab] = ^/'^SnABNqq, and Nqq is a 
zero mode. P^ has eigenvalues ±1, and is diagonalizable 
for any d. 

Putting together these results, we see that the u evolu- 
tion system is completely ill-posed for d = (Garfinkle's 
version of the equations), but is strongly hyperbolic for 
d < —1/2. The most natural value is d — —2, for which 
all characteristic speeds are (0, ±1). 

In order to find the most general e that is conserved, 
we parameterize the most general scalar e and vector F" 
that are quadratic in u without the use of background 
fields (in particular must not appear). We then use 
dt€ ~ 9qP" to determine their coefficients. The result 
for d = — 2 is 



Co 



a/3 



N 



2 



pa 



+ Ciri + C2 [Ea') + C3 {N'^aY 

3co [^ap{r^ - AP) + eaM^^^^^] 



(42) 
(43) 



Clearly this is positive definite for cq, ci, C2 > and C3 > 
0. 

We have therefore shown that with the modification 
d = —2, Garfinkle's version of Uggla et al's equations 
is mixed symmetric hyperbolic/parabolic, and from the 
theorem we have cited the resulting initial value problem 
is well-posed. 



IV. CONSTRAINT PROPAGATION 

We should stress that what we have proved is that 
the Cauchy problem for the evolution equations ("free 
evolution") is well-posed independently of whether the 
initial data obey the constraints or not. This is the well- 
posedness result required for stability of numerical free 
evolution schemes, because numerical error will always 
generate small constraint violations even if the initial 
data obey the constraints. 

Nevertheless, the constraints are compatible with the 
evolution equations in the sense that if we denote the 
constraints schematically by a vector c, the expression 
for dtc obtained by substituting the evolution equations 
is homogeneous in c and daC, so that c{x, 0) = im- 
plies dtc{x, 0) = 0. In order to guarantee a well-posed 
continuum problem solving all Einstein equations (evolu- 
tion and constraints) we should check that the constraint 
evolution is also well-posed. Together with its closure 
this guarantees that c(x, t) = is the unique solution of 
c(x,0) = 0. 

It is true in general that if the main system of evolution 
equations is strongly hyperbolic and the constraint evolu- 
tion system closes, then the constraint evolution system 
is also strongly hyperbolic .4]. However, there is no such 
theorem for a mixed parabolic/hyperbolic problem, and 
therefore we carry out the analysis explicitly. The prin- 
cipal part of the constraints is 



{Ccom)ap 


— 29[Qi?^], 


(44) 


Cg 


^ \daA'^, 


(45) 


{CcT 




(46) 


Cq 


c q+^dar^, 


(47) 


{CjT 


~ dpN'^'^ + e"'^''df!A^, 


(48) 


[CwT 




(49) 



One may wonder at the appearance of the undifferen- 
tiated q in the principal part of Cq, but its presence 
is required if the principal part of the constraints is to 
close under the principal part of the evolution equations. 
These are given by and l(77|) - (|^ . with the excep- 
tion that including the cross terms in the principal part 
we must replace by 



dtVa ^ daq- 



(50) 
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We find that under the principal part of the evohition 
equations the principal part of the constraints evolves as 



at(Ccom)l/3 - 0, (51) 

dtCo ^ I (rf+0a4Cc)", (52) 

dtiCcr ^ ie^/'^a^ [(Cj)^ - (Cw)^] 

- 9"(CG-2Cq), (53) 

dtC, ~ 0, (54) 

dtiCjr ^ de^^-'dpiCc)^, (55) 

dt{CwT - 0. (56) 



An analysis similar to the one carried out for the main 
system shows that the constraint system is strongly hy- 
perbolic for d < —1/2, with the same speeds that arise 
in the main system. The constraint evolution system is 
therefore well-posed if the main system is. With hind- 
sight the fact that the constraint system does not have 
a parabolic part can be explained by the fact that the 
parabolic part of the main system is the evolution of the 
slicing, that is, gauge, which should not affect the con- 
straint evolution directly. 

V. CONCLUSIONS 

We now consider the implications of this result for nu- 
merical simulations of the systems discussed. For an ill- 
posed system, unbounded rates of growth are achieved in 
the limit as wavelength goes to zero. This is what occurs 
in the scalar sector of our system for d > —1/2. 

However, the Einstein equations are well-posed physi- 
cally in the sense that they are well-posed in appropriate 
variables, and so the instability can only be due to a 
gauge transformation or to constraint violation. In the 
high frequency limit the constraints (|17I19I) imply that 
daS" = and da'S,"^ = 0. Therefore, if Ua is chosen 
to be the direction of the wave number of the mode, the 
components s„ and (as well as Syi„) vanish if the 
constraints vanish. Therefore all scalar modes, and hence 
all unstable modes for d > —1/2, are constraint violating 
modes. 



In a numerical simulation, the effective wavelengths 
involved are not arbitrarily small; but are instead lim- 
ited to a size of order the spatial grid spacing. The rate 
of growth of the constraint violating mode will then be 
limited by the spatial resolution. Furthermore, since the 
constraints will be enforced (up to numerical truncation 
error) in the initial data, the initial amplitude of the un- 
stable mode will be small. 

In numerical simulations of the d — system (which 
we have shown here to be ill-posed) are performed with 
a numerical resolution of 50 grid points for each spatial 
dimension. We have performed numerical simulations 
of the same initial data on the same spatial grid for the 
d — —0.6 system (which is well-posed). The results of the 
two simulations are effectively identical. The reason for 
this is that the nonlinear terms in the evolution equations 
tend to damp constraint violations. Thus in the d = 
evolution the unstable constraint violating mode (which 
starts out at very low amplitude) does not have time to 
grow to an appreciable amplitude before it is damped 
by the nonlinear terms. (The reason for using d = —0.6 
rather than d = —2 is that in the d — —2 system lower- 
order terms can lead to the growth of constraint violat- 
ing modes). However, we have also performed numerical 
simulations of the the d = and d = —2 systems with a 
spatial resolution of 1200 gridpoints. We have done this 
by choosing initial data that depends on only one spa- 
tial coordinate. Here, the results of the two simulations 
are very different: the d = ~2 simulation yields stable 
evolution, while the d = simulation has an unstable 
constraint violating mode that grows from a very small 
amplitude to an appreciable one. Thus, for sufficiently 
good resolution, the modification of our work to the evo- 
lution equations of is needed for a stable numerical 
evolution. 
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